Integration

Introduction

In the back of your textbook, there are many integration formulas. In times long past, before mathematics software existed, there were entire multivolume books of these formulas, and mathematicians had to either use these books or derive everything by hand. Mathematics software makes the problem of computing integrals much, much simpler, as all of these formulas are made available to the software - you don't have to go and leaf through a collection of tables to find the answer.

Our goal with this notebook is to understand how Mathematica approaches the problem of integration. We will use Mathematica to compute symbolic integrals, and will learn to generalize these results into integration formulas of our own. We will also examine some simple methods for computing integrals numerically.

Symbolic Integration

Mathematica's Integrate command is very versatile. If we give it an expression and a variable to integrate by, it will give us the indefinite integral. Note that Mathematica does not include the constants of integration.

Integrate[Log[x],x]

Alternatively, we can give it a range in order to compute definite integrals:

Integrate[Log[x],{x,1,5}]

We can even use Integrate as part of another expression:

Solve[{Integrate[x^3-2x+1,{x,1,p}]==5},p]

Sometimes, Integrate makes assumptions that we are not making. Consider the following integral:

Integrate[ArcSin[x],{x,Pi,2*Pi}]

A Calculus student is expected to give a very different answer from the output of the above command.

That being said, we can still make Integrate work for us. One of the ways of doing this is by looking at the output of simple formulas and looking for patterns. We saw an indefinite integral of ln(x), what about x ln(x)?

Integrate[x*Log[x],x]

It should look somewhat similar to the indefinite integral of ln(x). Let's try some higher powers:

Integrate[(x^2)*Log[x],x]
Integrate[(x^3)*Log[x],x]
Integrate[(x^4)*Log[x],x]

In general, it seems like we have a formula that looks like this: ∫(x^n)ln(x)dx=(((x^(n+1))ln(x))/(n+1))-((x^(n+1))/(((n+1))^2)).

Integrate[(x^n)*Log[x],x]
Simplify[Integrate[(x^n)*Log[x],x]==(x^(n+1)*Log[x]/(n+1)-x^(n+1)/(n+1)^2)]

This is an example of one of the most powerful formulas, the Integration By Parts formula ∫udv=u v-∫vdu. Many, many other formulas are simply applications of Integration by Parts.

Another important technique is called Partial Fractions Expansion. Let's look at this integral:

Integrate[1/(x^2-1),x]

The answer is a sum of two terms. If we try to go back a few steps into how to get this answer, then that means at some point we should be able to break up the expression (1/((x^2)-1)) into two parts. What's more, the fact that the terms are logarithms means that these terms prior to integration should look something like (1/x). The purpose of Partial Fractions Expansion is to break up large polynomials that are in the denominator of a rational expression. Mathematica has a function called Apart that can perform this expansion for us:

Apart[1/(x^2-1)]

Integrating the individual terms does indeed give us the antiderivative that we found previously.

Apart also works on more complicated expressions:

Apart[(2*s+2)/((s^2+1)*(s-1)^3)]

Numeric Integration

Sometimes Integrate doesn't give us useful output. It's just not smart enough to solve the integral for us:

Clear[x]
Integrate[(1+Log[x])*Sqrt[1+(x*Log[x])^2],{x,1,10}]

In times like these, we have to approximate the value of the integral numerically.

Before we can understand how to approximate the value of a definite integral with Mathematica, we need to first understand some concepts from computer science:

Conditional Statements

A conditional statement is a way to express decision making. We use one expression which can take the value "True" or "False," and based on that expression, we evaluate other expressions. In Mathematica, the simplest conditional statement is the If function. If takes three arguments: the first argument is the expression to decide with, the second argument is the input to evaluate if the expression is true, and the third argument is the input to evalaute if the expression is false.

1>2 
2>1
If[1>2,2+2,2+3]
If[2>1,2+2,2+3]

Another important conditional statement is Piecewise. It allows us to construct piecewise-defined functions. The way to use Piecewise is to give it a list, and each element of the list is an expression-condition pair.

f[x_]:=Piecewise[{{x,x<1},{x^2,x≥ 1}}]
Plot[f[x],{x,-1,2}]

The most useful expressions that can take values of "True" or "False" available to us are called relational operators. They allow us to test the relationship between two other terms. We've seen > (greater than), < (less than), and >= (greater than or equal to) already. Some other important ones:

If[1==2,2+2,2+3] (*Equal To (==)*)
If[1≠ 2,2+2,2+3] (*Not Equal To (!=)*)
If[1≤ 2,2+2,2+3] (*Less Than Or Equal To (<=)*)

Looping

If we want to evaluate many very similar expressions, we should use looping to accomplish the task. A loop consists of Mathematica commands to evalaute "inside" of the loop, as well as a condition which is used to decide when the loop ends. The easiest type of loop is the Do loop. It takes an expression and a variable range, and evaluates the expression once for each integer in the range:

Do[x=i;Print[x],{i,1,10}]

A more flexible kind of loop is the While loop. It is given a condition first and a command or sequence of commands second. It evaluates the commands until the condition is "False."

x=1
While[x<5,x=x+1;Print[x]]

Be careful: it is possible to make an infinite loop with While. If this happens, you need to quit the kernel. The command to do this is found in the Evaluation menu.

x=1
While[x<5,x=x-1;Print[x]] (*Quit the kernel after running this!*)

In order to make our conditional statements and loops easier to read, it is useful to spill out a single command over multiple lines:

Do[
x=i;Print[x],
{i,1,10}
]

Our goal is to approximate the value of an integral (∫_a^b)f(x)dx, for some function f(x) and interval (a,b). The simplest methods divide up this interval into n subintervals, each of length h=((b-a)/n).

Left and Right Hand Rules

One way is to use the fact that an integral is defined as a Riemann sum. Page 532 of your textbook describes the left-hand rule:

(∫_a^b)f(x)dx≈(∑_(i=1)^n)f((Subscript[x, i-1]))h

For the function f(x)=(e^x)+3x, a=1, b=3, and n=10, we can implement this with loops as follows:

Clear[x]
f[x_]:=Exp[x]+3*x
a=1;b=3;n=10;h=(b-a)/n;sum=0;
Do[sum=sum+f[a+(i-1)*h]*h,{i,1,n}]
N[sum]
N[Integrate[f[x],{x,a,b}]]

Notice how the value we get is only approximate. The first number is what we get from the left hand rule. The second number is what we get from evaluating the integral symbolically first, and then converting that to a decimal answer.  We've seen the power of N before, so we can expect the second number to be very close to the true answer.

We can also compute the left-hand rule as a sum using the Sum command:

Sum[f[a+(i-1)*h]*h,{i,1,n}]
N[%]

We can use the right-hand rule to get similar results:

(∫_a^b)f(x)dx≈(∑_(i=1)^n)f((Subscript[x, i]))h

Trapezoidal Rule

A slightly better method is to average the left-hand and right-hand rules to get the trapezoidal rule, described on page 533 of your textbook:

(∫_a^b)f(x)dx≈((f(a)+f(b))/2)h+(∑_(i=1)^(n-1))f((Subscript[x, i]))h

Computing this with loops is roughly the same complexity:

Clear[x]
f[x_]:=Exp[x]+3*x
a=1;b=3;n=10;h=(b-a)/n;
sum=(f[a]+f[b])*h/2;
Do[sum=sum+f[a+i*h]*h,{i,1,n-1}]
N[sum]
N[Integrate[f[x],{x,a,b}]]

Notice how the value we get with the trapezoidal rule is still off, but is much closer than what we get through the left- or right- hand rules.

We have to add the extra terms when performing a sum, as well:

Sum[f[a+i*h]*h,{i,1,n-1}]+(f[a]+f[b])*h/2
N[%]

Simpson's Rule

We can get an even better approximation by using Simpson's rule, which is given on page 538 of your textbook:

(∫_a^b)f(x)dx≈((f(a)+f(b))/3)h+(∑_(i=1)^(n/2))f((Subscript[x, 2i-1]))((4h)/3)+(∑_(i=1)^(n/2-1))f((Subscript[x, 2i]))((2h)/3)

Implementing this with loops looks like the following:

Clear[x]
f[x_]:=Exp[x]+3*x
a=1;b=3;n=10;h=(b-a)/n;
sum=(f[a]+f[b])*h/3;
Do[sum=sum+f[a+(2*i-1)*h]*(4*h/3),{i,1,n/2}]
Do[sum=sum+f[a+2*i*h]*(2*h/3),{i,1,n/2-1}]
N[sum]
N[Integrate[f[x],{x,a,b}]]

Implementing Simpson's rule with sums looks like this:

(f[a]+f[b])*h/3+Sum[f[a+(2*i-1)*h]*(4*h/3),{i,1,n/2}]+Sum[f[a+2*i*h]*(2*h/3),{i,1,n/2-1}]
N[%]

Monte Carlo

The Monte Carlo method is not typically used with the types of problems you will encounter in Calculus classes, but it does have applications when the region of integration is strange. The idea is that the average value of our function over the interval can be expressed as (((∫_a^b)f(x)dx)/(b-a)). The details of this calculation can be found on page 375 of your textbook. The steps to perform in order to approximate the integral are the following:

Generate n random values in the interval [a,b].

Evaluate f at these random values.

Average these values and then multiply by b-a.

For the same example we've been using thus far, we can implement this as follows:

Clear[x]
f[x_]:=Exp[x]+3*x
a=1;b=3;n=10;sum=0;
Do[random=RandomReal[{a,b}];
sum=sum + f[random],
{i,0,n-1}]
sum=sum*(b-a)/n;
N[sum]
N[Integrate[f[x],{x,a,b}]]

An interesting side effect is that, since random numbers are being used, the value we get is also random. It could be very close to the exact value, or it could be quite far. Note that the answer is different every time you execute the above code.

NIntegrate

Mathematica has also has a built in command called NIntegrate that allows us to perform numeric integration. It is used just like Integrate, except that you must always compute definite integrals with it. Mathematica uses various different methods - most of them more complicated and more accurate than the simple methods explained above - to arrive at its answer. You might study such methods in a course on computational or numerical mathematics.

NIntegrate[f[x],{x,a,b}]

Exercises

Assignment 1: Understanding Antiderivatives

Use Mathematica to compute antiderivatives for the following functions:

y(x)=x sin x

y(x)=((arctan x)/((x^2)+1))

You should give at least two antiderivatives for each function.

Assignment 2: Working with Symbolic Integrals

Previously, we saw that attempting to compute ∫(1+ln(x))(Sqrt[1+(((x ln(x)))^2)])dx with Integrate didn't work. However, if we help Mathematica to get started, it may be able to finish off the rest of the work. Express the integral in a different form and use Mathematica to compute its antiderivative. Your answer should include:

A simplification of ∫(1+ln(x))(Sqrt[1+(((x ln(x)))^2)])dx  that works with Integrate.

An antiderivative of  (1+ln(x))(Sqrt[1+(((x ln(x)))^2)]).

Assignment 3:  Obtaining an Integration Formula

Use Mathematica to evaluate the first few integrals of the form ∫(1/(x^n))ln(x)dx, n≥ 2 and generalize this result into an integration formula. Your answer should include:

Antiderivatives for ∫(1/(x^2))ln(x)dx,∫(1/(x^3))ln(x)dx, and ∫(1/(x^4))ln(x)dx.

The general formula for ∫(1/(x^n))ln(x)dx, n≥2.

A comparison between your formula for ∫(1/(x^n))ln(x)dx and Mathematica's output for this same integral.

Assignment 4:  The Effect of the Number of Points

Investigate the effect of the number of subintervals n in the Left-Hand rule, the Trapezoidal Rule, and Simpson's rule on the quantity obtained by the difference between the result of the algorithm and the true value of the integral. Do this for the integrals (∫_0^(Π/2))((sin(x))/(1+(((cos(x)))^2)))dx and (∫_0^5)(x^23)(e^x)dx. Your answer should include:

The differences between the approximate values of the integrals for n=10,20,30, and the true value of the integrals. (6 answers for each method listed above)

A statement of how fast each method improves as n increases.